Simulating Impacts of Precipitation on Ice Cover and Surface Water Temperature Across Large Lakes

2.1 Ice-Hydrodynamic Model

We use the unstructured grid, Finite Volume Community Ocean Model (FVCOM, Chen et al., 2006, 2013) to simulate the Great Lakes hydrodynamics. FVCOM is a three-dimensional, free-surface, primitive equation, sigma-coordinate oceanographic model that solves the integral form of the governing equations. FVCOM has been applied in several studies of the coastal ocean, including successful application to operational forecasting in the Great Lakes (E. J. Anderson et al., 2010, 2015, 2018; E. J. Anderson & Schwab, 2012, 2013; Bai et al., 2013; Niu et al., 2015; Luo et al., 2012). In this work, the model is configured separately for Lake Superior, Lake Erie, and Lake Ontario, while Lake Michigan and Lake Huron are handled by the same model as the lakes are connected by the Straits of Mackinac and form a single system. Horizontal grid resolution in each model ranges from roughly 200 m near the shoreline to 2,500 m offshore, with 21 vertical sigma layers evenly distributed throughout the water column. As a result, the numbers of triangular elements in the models are roughly 23,000 for Lake Superior, 31,000 for Lake Michigan-Huron, 12,000 for Lake Erie, and 35,000 for Lake Ontario. For Lake Erie, the grid configuration is based on the National Oceanic and Atmospheric Administration (NOAA)'s Lake Erie Operational Forecast System (Anderson et al., 2018), and the remaining lakes have similar resolution and configuration. Horizontal diffusion is handled by the Smagorinsky parameterization (Smagorinsky, 1963), and vertical diffusion is handled by Mellor-Yamada Level-2.5 turbulence closure scheme (Mellor & Blumberg, 2004; Mellor & Yamada, 1982) with surface wave breaking parameterization by Craig and Banner (1994). Turbulent latent and sensible heat fluxes (not associated with precipitation) are calculated from the Coupled Ocean-Atmosphere Response Experiment (COARE, Fairall et al., 1996, Fairall et al., 1996, Fairall et al., 2003) algorithm. Modeled depths (Figure 1) are interpolated from 3-arc sec bathymetry data from the NOAA National Centers for Environmental Information.

An unstructured grid version of the Los Alamos Sea Ice model (CICE; Gao et al., 2011; Hunke et al., 2015) has been coupled within FVCOM. The CICE model includes components for ice thermodynamics and ice dynamics, using elastic-viscous-plastic rheology for internal stress (E. C. Hunke & Dukowicz, 1997) and produces two-dimensional fields of ice concentration, thickness, and velocity. A multicategory ice thickness distribution model (Thorndike et al., 1975) is employed in CICE to represent the sub-grid-scale distribution of ice thickness in response to mechanical and thermal forcing. In this study, five categories of ice thickness are defined (5, 25, 65, 125, and 205 cm). The ice surface albedo depends on surface temperature, snow depth, and thickness of ice, as well as the visible and infrared spectral bands of the incoming solar radiation (Briegleb, 1992). Specifically, visible and near-infrared spectral albedos are distinguished for snow and ice. For the bare ice (snow) albedo, 0.78 (0.98) and 0.36 (0.70) are used for visible (<700 nm) and near-infrared (>700 nm) wavelengths, respectively. The ice (snow) albedos in both spectral bands decrease by 0.075 (0.15) as the surface temperature rises from −1°C to 0°C. As ice thickness decreases from 50 cm to 0, the ice albedo decreases smoothly (using an arctangent function) to the open water albedo, 0.06. The total albedo is an area-weighted average of the ice and snow albedos, where the fractional snow-covered area f_snow is defined as f_snow = h_s/(h_s + h_snowpatch), where h_s is snow depth and h_snowpatch is a constant set at 2 cm. The snow and ice albedo parameterization is consistent with the one in the previous version of the Community Climate System Model (CCSM3) and is detailed in Hunke et al. (2019). For the density of snow ρ_s, a constant value of 330 kg/m³ was used. The thermal conductivity of snow and ice were set to 2.03 and 0.30 W/m/K, respectively.

At ice-covered cells, the net momentum transfer is calculated as a weighted average of the air-water and ice-water stresses by areal fraction of ice. The air-ice drag coefficient C_{D_ai} is a function of wind speed U, given as C_{D_ai} = (1.43 + 0.052 U) × 10⁻³ and the ice-water drag coefficient is 5.5 × 10⁻³. Similarly, the net heat transfer is calculated as a weighted average of the air-water and ice-water heat fluxes. The ice-water heat fluxes are calculated based on the bulk transfer formula (Maykut & McPhee, 1995). The formation of snow ice, which forms when flooded snow refreezes, is parameterized by calculating the depth of snow submerged below the water surface based on the hydrostatic assumption (CICE; Consortium, 2019). The submerged snow is converted to ice.

For this particular study, a few updates were made on the FVCOM + CICE codes: First, the parameterization for the air-water drag coefficient was updated. In the original FVCOM code, the air-water drag coefficient is calculated as a function of wind speed (Large & Pond, 1981) and does not vary with stratification in the surface boundary layer. In this study, we used the wind stress (momentum flux) calculated by COARE (within FVCOM+CICE) to incorporate the dependency on atmospheric stratification. Second, calculation of the heat fluxes due to precipitation was added. The sensible heat flux from precipitation was implemented using the existing method in COARE (Gosnell et al., 1995). As for the latent heat flux from precipitation (cooling due to snow melting), it was calculated as the linear function of snow precipitation (Duffy & Bennartz, 2018). In these parameterizations, the sensible and latent heat flux due to precipitation, H_sp and H_lp are described as bulk formulae:

(1)

(2)

where ρ_w is the density of water and c_pw and c_pi are the specific heats of water and ice, respectively. L_i is the latent heat of melting for ice, and P is the rate of precipitation. T_sfc and T_precip are the water surface temperatures and the temperature of rain droplet or snow flake, which is approximated as the wet-bulb temperature T_b. The wet-bulb temperature was derived based on the psychrometer equation and its derivation is based on Gosnell et al. (1995). In equation 1, H_sp can be both for rain and snow and in equation 2, H_lp is for snow only. Precipitation was considered as snow when the wet-bulb temperature T_b is below −2°C and otherwise rain. For H_sp and H_lp, the sign conventions are downward positive (i.e., warming of lake surface) and upward negative (i.e., cooling of lake surface). H_sp and H_lp were only for the water surface. If a model cell is covered with lake ice, H_sp and H_lp were scaled with a fraction of open water. Note that ρ_w instead of ρ_i (the density of ice) is used in equation 2 as the precipitation rate P is based on the liquid phase of water.

The selected experiment period was from 11 November 2014 to 1 November 2016. This period covers one severe winter with high ice cover (2014–2015) and a mild winter with low ice cover (2015–2016) and is ideally suited to evaluate the impacts of precipitation on variable ice conditions, as well as thermal conditions afterwards. The control case started on 1 April 2014 with uniform 4°C temperature and static conditions (cold start option) and forced by interpolated observations from coastal and offshore meteorological stations that are corrected for overwater conditions (Beletsky et al., 2003) until September 2014, after which an operational weather forecast from the High-Resolution Rapid Refresh (HRRR, described in section 2.2) was utilized as the forcing data. All the three experiments, which are described in section 2.3, were started from a restart file of this control case on 11 November 2014. This start date for the three numerical experiments allowed to cover the first winter storm over the Great Lakes region in the winter of 2014–2015 (Fujisaki-Manome et al., 2017) when notable snow and therefore was appropriate to examine the impacts of precipitation (snow or rain) on lake ice and water temperature.

2.2 Data

For the meteorological forcing to the FVCOM-CICE model, the HRRR (Benjamin et al., 2016,Benjamin et al., 2016) was selected. HRRR is a relatively new, cloud-resolving, convection-allowing weather model with horizontal resolution of 3 km that became operational in September 2014. Radar data are assimilated in the HRRR every 15 min over a 1-hr period, adding further details to those provided by the hourly data assimilation from the 13-km radar-enhanced Rapid Refresh. The meteorological forcings from HRRR at forecast hour 2 were applied to the FVCOM-CICE model on an hourly basis by providing wind speed at 7.8-m height, surface air temperature at 2-m height, relative humidity at 2-m height, and precipitation rate. The COARE algorithm used in this study is able to take account of these reference heights to derive the momentum and turbulent heat fluxes in the surface boundary layer.

Cross-evaluation of precipitation rate from the HRRR was conducted beforehand in comparison with the merged product of the Canadian Precipitation Analysis (CaPA) and the Multisensor Precipitation Estimate (MPE). CaPA is a real-time gridded precipitation product provided by Environment and Climate Change Canada and covers all of North America. MPE is the composite analysis of weather radars, rain gage data, and satellite precipitation estimates and is provided by the U.S. National Weather Service. The merged CaPA-MPE product combines these two data sets into 10-km gridded data and has spatially continuous precipitation estimate over the Great Lakes. The merged CaPA-MPE product is widely used by the binational water balance community for hydrologic analysis and lake management and therefore an adequate data set to compare with HRRR's overlake precipitation to ensure its validity. Further details can be referred to Gronewold et al. (2018). We calculated daily overlake precipitation for each lake from the HRRR and compared with those from the merged CaPA-MPE product (Figure 2). The two data sets generally agreed well. Slightly larger values in HRRR, represented by a slope value s > 1, are likely because the CaPA-MPE rely heavily on radar estimates from the coastlines, which are sometimes limited in terms of coverage and obstacles (e.g., the radar at Marquette, Michigan, subject to the blockage by the Keweenaw Peninsula), and therefore could miss some overwater precipitation. Overall, the reasonable agreement of HRRR's overlake precipitation estimate with the merged CaPA-MPE provides confidence in using overlake precipitation from HRRR. It should be reminded that precipitation from HRRR has other advantages such as high spatiotemporal resolution and physically consistent representation based on the single model.

For the evaluation of water surface temperature and ice coverage, the Great Lakes Surface Environmental Analysis (GLSEA; https://coastwatch.glerl.noaa.gov/glsea/doc/) was used along with ice concentration maps from the National Ice Center. GLSEA is the composite analysis from six satellite measurements based on the Advanced Very High Resolution Radiometer and the Visible Infrared Imaging Radiometer Suite. The water surface temperatures are updated daily using information from the cloud-free portions of the previous day's satellite imagery. For pixels where data are not available, an interpolation algorithm is applied using cloud-free data over a ±10-day window. Therefore, the accuracy tends to be lower on cloudy days.

2.3 Numerical Experiments

We conducted three numerical experiments to assess the impacts of precipitation on lake ice and water temperature. In the control experiment (Expt. 1), no precipitation is considered. In the first precipitation experiment (Expt. 2), precipitation rate from the HRRR was applied to the model as the input data. The sensible and latent heat fluxes from precipitation (H_sp and H_lp) were calculated by the models in these experiments. The other heat flux components (i.e., shortwave radiation, longwave radiation, turbulent sensible, and latent heat fluxes) were calculated dynamically in the model and H_sp and H_lp are aggregated with the other heat flux components into the net heat flux. We also conducted the supplemental precipitation experiment (Expt. 3), where the other heat flux components than H_sp and H_lp were prescribed from the outputs of Expt. 1. The net heat flux over water was extracted from the Exp. 1 (no precipitation) results, and this net heat flux was fed to Expt. 3, where the model calculated the heat fluxes only for H_sp and H_lp and added them to the prescribed net heat flux. The purpose of Expt. 3 was to evaluate the maximum potential of H_sp and H_lp impacts on the water temperature. In Expt. 2, both water temperature and the other heat flux components were allowed to respond to the precipitation heat fluxes (i.e., perturbation), while in Expt. 3, only water temperature was allowed to respond to the perturbation because the other heat flux components were fixed. Therefore, the impacts on water temperature are meant to be maximized in Expt. 3.

3.1 Verification of Modeled Ice Conditions

Seasonal evolutions of ice coverage were reasonably simulated for each of the Great Lakes by the models (Figure 3) in comparison with the National Ice Center (NIC) analysis. This is consistent with the long-term model validation for Lake Erie and Lake Michigan-Huron by Anderson et al. (2018). The contrasts of high and low maximum ice coverage reflected the anomalously cold (2014–2015) and warm (2015–2016) winters. The spatial patterns of ice concentration at the early, middle, and late ice stages agreed with the NIC analysis (Figure 4). There were some overestimations of ice coverage during late ice season (March–April) in 2015, which were most notable in Lakes Superior and Huron and in Lakes Erie and Ontario (Figure 3). With precipitation (Expts. 2 and 3), ice coverage declined slightly faster than that in Expt. 1. The root-mean-square errors (RMSEs) of ice coverage were slightly decreased in Expt. 2 from the results in Expt. 1 during the winter of 2014–2015, except for Lake Ontario (Table 1). In the low-ice winter of 2015–2016, the RMSEs were similar among the three experiments. As anticipated, the simulated snow and ice conditions in Expt. 3 were nearly identical to that in Expt. 2. Therefore, in the following discussion in this section, we will focus on the results from Expts. 1 and 2 only.

Figure 5 shows the time series of simulated ice volume for each of the Great Lakes. A notable reduction in ice volume occurred by inclusion of precipitation (i.e., snow). The maximum ice thickness, which was obtained by dividing the ice volumes by the ice areas, was reduced by including precipitation (i.e., Expt. 2 versus Expt. 1) in all lakes for 2014–2015: 7, 6, 2, 8, and 1 cm for Lakes Superior, Huron, Michigan, Erie, and Ontario, respectively. In the low-ice winter of 2015–2016, the reduction in ice volume was less evident, but the maximum ice thickness was reduced by 3 and 1 cm in Lake Huron and Ontario, respectively. In the precipitation experiments, the shorter ice durations by 1–4 days were seen for most of the lakes (Table 2).

The modeled snow cover on the ice in the precipitation experiment (Expt. 2) was 0–8 cm in the winter of 2014–2015 and less than 1 cm in the winter of 2015–2016 (Figure 6). The observations of snow cover on the ice that can be compared with the model results are rare. There have been limited visual observations of snow depth on the ice in Lake Erie reported by the U.S. Coast Guard during their helicopter flyovers, in collaboration with NOAA Great Lakes Environmental Research Laboratory. These observations occurred on 1 or 2 days per year in late February or early March from 2008 to 2015 (with a few skipped years). In the winter of 2014–2015, the reported values on 25 and 26 February 2015 over five stations in Lake Erie ranged from 0 to 10 cm. The range does not conflict with the modeled values (Figure 6), but care should be taken as the range from the reports is for point values, while Figure 6 shows the lake-wide means. The ice surface albedo ranged from 6% (bare, very thin ice) to 60% in the control experiment. For comparisons, Bolsenga (1969) reported 10–46% for snow-free ice in Lake Superior, being consistent with the model results. The modeled albedo increased by 0–20% when precipitation was considered. This is consistent with Bolsenga (1987) who reported 62–80% for snow albedo in Michigan locations. By late March, snow melted away and the ice surface albedo went down quickly to be similar to that in Expt. 1 (Figure 6). The spatial patterns of snow depth presented high spatial heterogeneity (Figure 7), which also contributed to the spatial variation in the ice surface albedo.

3.2 Snow Impacts on Ice Conditions

Inclusion of snow reduced the total volume of ice despite that snow cover on the ice increased surface albedo in early winter of 2014–2015 (Figure 6), which reflected more sunlight back to the atmosphere and therefore prevented ice surface from melting. Note that snow also served as an input to the ice mass by formation of snow-ice. Apparently, the thermal insulation of snow cover on the ice (Sturm et al., 2002; Warren et al., 1999) overcame the impacts of increased surface albedo and snow-ice formation. In the precipitation experiment (Expt. 2), the surface albedo quickly reduced in March, when the atmosphere started to warm, and snow cover started to melt away. In this period, the surface albedo in the precipitation experiment was slightly below that in the control experiment (e.g., Lake Erie in the winter of 2014–2015, Figure 6). This was because the modeled lake ice was thinner in Expt. 2 than in Expt. 1. As described in section 2.1, the parametrized ice albedo reduced as ice became thinner (Briegleb, 1992). By including snow, the modeled ice thickness became thinner due to the dominating thermal insulation effect as described above. This resulted in the lower surface albedo once the ice lost snow cover on it. The shortened ice duration by inclusion of precipitation was likely because the reduced total ice volume by snow cover on the ice (manifestation of precipitation) allowed earlier delay of lake ice. Overall, inclusion of precipitation improved the models in reproducing ice areas and durations, as demonstrated by the reduction in RMSE of ice coverage.

In the version of model used in the study, the amount of formed snow ice was not directly outputted, not allowing the direct assessment. However, estimations of maximum possible lake-wide mean snow depth were made using the modeled freezing precipitation amount, that is, P in equation 2. The estimates resulted in 6, 3, 3, 2, 5 cm (4, 2, 5, 1, 2 cm) for Lake Superior, Lake Huron, Lake Michigan, Lake Erie, and Lake Ontario for the winter of 2014–2015 (2015–2016). If they were all to form snow ice, the estimated maximum possible additions to the total ice volume would be 1.7, 0.7, 0.7, 0.2, 0.3 km³ (1.2, 0.4, 1.1, 0.1, 0.1 km³). These numbers were fairly small compared with the range in Figure 5. Thus, inputs from snow ice formation were likely minor to the total ice volume. This is somewhat contradictory to the findings in the polar oceans (e.g., Fichefet & Morales Maqueda, 1999) and other lakes (e.g., Ohata et al., 2017), but in these study areas, snowfall amounts over the ice surface were much larger than the estimates in this study, providing favorable conditions for snow ice formation.

For comparison with classical freezing degree day (FDD) ice growth models, Tables 3 and 4 show the estimated maximum ice thickness by the equations proposed by Lebedev (1938), referred in Bilello, 1961) and Ashton (1989). Note that the Lebedev model was based on the sea ice observations in the Arctic shelf seas. The Ashton model was for lake ice but no snow cover was considered. Overall, in the winter of 2014–2015, the estimates by the simple FDD ice growth models were in agreement with the three-dimensional ice-hydrodynamic model results. The FDD ice growth models tended to provide higher ice thickness. This was pronounced in Lake Superior, Lake Huron, and Lake Michigan in the winter of 2015–2016 and in all lakes in the winter of 2016–2017. One obvious reason was heat (i.e., warm water) carried from the preceding summer, which prevented ice from growing in reality. An FDD ice growth model typically takes air temperature as the single input for its simplicity and therefore is limited in application to dimictic lakes like the Great Lakes. Snow inclusion in the Lebedev model resulted in some reductions of ice thickness estimates for a constant snow depth h_s of 6 cm, as was seen in the three-dimensional model results. However, with h_s = 3 cm, ice thickness estimates resulted in unrealistic increase from the no snow results in both winters. Again, the Lebenov model was an empirical model and the data fitted to the equation had h_s and FDD mostly greater than 10 cm and 1,000°C/day, respectively (Figure 11 in Bilello, 1961). Therefore, its application to the Great Lakes should be limited as h_s and FDDs are typically below these ranges in a normal winter. The Ashton model provided closer results to the three-dimensional model results than the Lebenov model did, but the higher ice thickness estimates were also pronounced in the winter of 2015–2016, likely due to not accounting of warm water carried from the preceding summer. Thus, the simple FDD ice growth models were useful for estimating ice thickness at the first order but tended to provide higher estimates than the three-dimensional model did due to lack of the thermal history in the lakes. It should also be noted that the three-dimensional model also takes account of mechanical thickening of ice, such as ridging and rafting.

3.3 Precipitation Impacts on Heat Fluxes and Water Temperature

The lake-wide mean precipitation heat fluxes H_lp and H_sp were in general intermittent due to their association with snow or rain events but presented expected seasonality (Figure 8). Snowfall during winter caused negative H_lp, (i.e., cooling of water surface). The daily values of H_lp were below −50 W/m² even during major snowstorm events. This range was smaller than the net heat flux in the similar timeframe (−500–0 W/m², Figure 9), however, not negligible. With a simple column equation, the water surface with a mixed layer depth of 10 m would be cooled at 0.1°C/day under a constant cooling of −50 W/m². Depending on the winter mixed layer depth, which can be shallower (deeper) than 10 m at low (high) wind conditions, the rate could be faster or slower. The sensible heat flux due to precipitation, H_sp, was overall smaller in terms of magnitude than H_lp. Positive H_sp was due to warm rain relative to the water surface from spring to summer.

Water surface temperature was reasonably simulated by the models in comparison with the analysis from GLSEA (Figure 10a). The contrast between the anomalously cold and warm winters was demonstrated in the minimum values of lake-wide mean water surface temperature: In the winter of 2014–2015, the minimum values were around the freezing temperature for all five lakes, while in the winter of 2015–2016, the minimum daily values ranged from 0.3°C (Lake Erie) to 2.3°C (Lake Ontario). Inclusion of precipitation resulted in only minor changes in lake-wide mean water surface temperature. Biases in the modeled lake-wide mean water surface temperature were very similar among the three experiments (Figure 10b). However, the three experiments presented minor differences in three occasions. First, although not clearly visible in the lake-wide mean time series in Figure 10b, major snowstorms caused intermittent, local cooling at the water surface. Figures 11a and 11b show an example of such cooling on 3 February 2015, after a major snowstorm over the Great Lakes region resulted in significant H_lp particularly over Lake Ontario and Lake Michigan. However, this signal was buried in the lake-wide average time series (Figure 10). Second, from late April to early May in 2015, there was slight warming in the two precipitation experiments relative to the control experiment (Expt. 1) in Lake Superior, Lake Huron, and Lake Erie. This coincided with the increase in the incoming shortwave radiation and earlier decay of ice cover in the precipitation experiments during the same time (Figure 9b and Table 2). The warming in Expts. 2 and 3 during this period was likely because the 1–4 days earlier ice-off allowed shortwave radiation to reach the water surface and therefore warmed the water temperature earlier. The warming was most notable in Lake Erie, Lake Superior, and part of Lake Huron, as in example snapshots on 16 April 2015 (Figures 11c and 11d). This impact was not directly due to H_lp or H_sp, but the indirect consequence of snow cover on the ice leading to the earlier decay of the ice. Third, from February 2015 to May 2016, the water surface temperature was slightly cooled in the precipitation experiments. This weak signal can be seen in the reduced biases for Lake Superior and Huron (Figure 10b). This coincides with the negative H_lp during the period (Figure 8a). This feature was amplified in Expt. 3, where the response to the precipitation perturbation was constrained to water temperature and therefore maximized. The cooling was most evident in Lake Superior, Lake Huron and Lake Michigan, as in example snapshots on 30 May 2016 (Figures 11e and 11f). The lake-wide mean water surface temperature could have been more sensitive to H_lp in the winter of 2015–2016 due to the low ice coverage during the season. The amplitudes of cooling in Figure 11 appear to be at the same order as the column model estimate by H_lp provided earlier in this study. However, even so, the H_lp impacts on the water surface temperature appear to be small compared with the biases against the observation and the seasonal variability, given that these signals are hard to see in the lake-wide mean time series in Figure 10.

The RMSEs of water surface temperature were improved in the precipitation experiments, but only slightly (Table 5). Thus, the overall precipitation impacts on the water surface temperature were found to be local and intermittent and be very minor from a lake-wide average perspective.

3.4 Model Uncertainties

Major uncertainties in this study were from overlake precipitation estimates and representation of snow cover on the ice. Even though we used the “best” available estimates of overlake precipitation, we consider that the observations included higher uncertainty than the other meteorology variables due to its extremely high variability in space and time. The threshold to distinguish snow and rain (i.e., the wet-bulb temperature −2°C, see section 2.1) was empirical (Knox et al., 2017; Olsen, 2003) although it should serve as a first-order accurate threshold at least. It should be noted that the current HRRR product provides snowfall. However, the authors did not have access to the variable in the archive to cover the entire simulation period.

Like in many other models, representation of snow cover on the ice accompanied various approximations in this study. These include snow blown by winds, aging of snow, and the associated changes in the density and heat conductivity of snow. The snow-ice formation parameterization is based on the hydrostatic balance and highly simplified. In reality, the flooding of snow-ice interface can occur more frequently in a wave environment, which could allow more formation of snow ice. Relative importance of these processes is still unclear, not only in the Great Lakes but also in the polar regions (Lei et al., 2016, 2017; Webster et al., 2018). In order to determine which processes are significant in controlling Great Lakes ice cover (and therefore a priority for modeling), more mature observations are needed to examine models. In particular, spatial observations for ice thickness, snow depth on the ice, and snow/ice albedo by satellite or aircraft missions would be valuable to aid in understanding lake ice mass, energy budget, and the roles of snow in these components. Process-oriented observations would also greatly help to inform the development of lake ice models, including both numerical models and simple FDD ice growth models.